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In a tight binding framework, we analyze the characteristics of electronic states in strongly dis- 
ordered materials (hopping sites are placed randomly with no local order) with tunneling matrix 
elements decaying exponentially in the atomic separation with various decay ranges I examined. 
We calculate the density of states (DOS) and the Inverse Participation Ratio (IPR) for amorphous 
atomic configurations in one, two, and three dimensions. With a finite size scaling analysis of the 
IPR statistical distributions, it is shown that states are either extended or localized for a particular 
energy, and phase portraits for wave functions are obtained showing extended and localized behavior 
in the thermodynamic limit. While we conclude that all states are localized in ID, in the 2D case 
there is a threshold for I above which some eigenstates appear to be extended and below which 
wave functions are entirely localized. For 3D geometries, there are two mobility boundaries flanking 
an intermediate range of energies where states are extended with eigenstates localized for energies 
above or below this range. While a zone of extended states persists even for very short I, the width 
of the region tends to zero exponentially (i.e. scaling as e _A// ') for very small decay length scales. 
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I. INTRODUCTION AND THEORETICAL 
FRAMEWORK 

A periodic crystal in the absence of disorder sup- 
ports extended Bloch waves within the bounds of energy 
bands 1 . However, the nature of the electronic states in 
a strongly disordered (i.e. amorphous) material where 
symmetry with respect to discrete translations is absent, 
is a more subtle question. Our aim is to examine the ef- 
fect of very strong disorder on transport characteristics. 

Operating in the framework of a tight binding model 
where electrons are localized in atomic sites or dopant 
impurities, we examine amorphous materials in one, 
two, and three dimensions where the locations of atoms 
are taken to be uncorrelated and randomly distributed 
within the medium. Atomic or dopant configurations 
with no local ordering of the sites are described as gas- 
like disorder— with relevance to the transport character- 
istics associated with expanded alkali metals^ as well as 
impurity bands in silicon. The characteristics of exciton 
states with respect to localization have been examined 
in the context of similar types of disorder™. In broader 
generality, formal analytical and computer studies have 
calculated the Density of States in amorphous materials 
with no correlations among the site positions 5-8 

Disorder, even in regular lattices, may be manifest 
as random site energies which can disrupt the extended 
character of itinerant states and thereby create condi- 
tions for localization. Our aim is to examine strongly 
disordered materials and the properties of the associated 
electronic states with respect to localization. However, 
we do not introduce a random site energy, and in this 
sense our work is complementary to studies where ran- 
dom potentials are superimposed on sites in a periodic 
crystalline geometry^—. Instead of examining a sys- 
tem on a regular lattice geometry, we calculate electronic 



wave functions for the fully amorphous case and exam- 
ine how strong positional disorder, in conjunction with 
tunneling matrix elements which decay exponentially in 
the separation of neighboring hopping sites, affects the 
characteristics of eigenstates with respect to localization. 

Off-diagonal disorder which enters in random varia- 
tions in tunneling matrix elements controlling hopping 
among sites in the absence of a random on-site po- 
tential, has been of interest since an early study by 
F. J. Dysoniir— . Calculations related to the density of 
states^— and aspects of localizatio n 4 ' 5 ' 21 have been car- 
ried out in the context of a 3D gas-like tight binding 
model. We report on results relevant to the aim of find- 
ing out by direct calculation in a large-scale statistical 
study the extent to which wave functions are localized in 
amorphous materials and under what conditions eigen- 
states are extended. Moreover, we obtain phase portraits 
showing for different energies and ranges / of the hopping 
integral domains of extended and localized states. 

In this work, we use for the tight binding Hamiltonian 
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where the sum over the index "i" ranges over the N parti- 
cles contained in the simulation volume, we take the hop- 
ping parameter to to be 3.0 electron volts, and the factor 
of "1/2" compensates for multiple counting of the hop- 
ping terms between atoms. The creation and destruction 
operators c' and c create and destroy occupied electronic 
orbitals at sites indicated by the subscript. For the hop- 
ping integral, we use V(rij) = e~ irij ^ s , where is the 
separation between sites i and j, s = p~ x l D is the typical 
inter-orbital separation, p is the volume density of sites, 
D is the dimensionality of the system, and 7 is a dimcn- 
sionless parameter. With the hopping integrals expressed 
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in this manner and the length scale for the decay of the 
hopping matrix element thus being I — s/7, large/small 
7 values correspond to decay lengths small/large in rela- 
tion to the typical interatomic separation. For the sake of 
convenience, we rescale coordinates such that p = 1, with 
a simple inverse relationship I = 7" 1 among the hopping 
length scale I and decay parameter 7 and V{rij) = e^ 71 "^ 
for the dependence of the tunneling matrix element. Al- 
though the hopping integral V(rij) is finite in range by 
virtue of the exponential decay, we nevertheless take into 
consideration hopping among all pairs of orbitals con- 
tained in the system of N sites with only a negligible 
increase (i.e. a contribution on the order of N 2 ) where 
direct diagonalization, which scales as A 3 and constitutes 
the main computational bottleneck, is used to calculate 
the electronic states. 

Our Hamiltonian does not incorporate perturbations 
of the on-site energy term (for convenience all site ener- 
gies are set to zero), and disorder instead enters in the 
off-diagonal tunneling matrix elements. Thus, our pro- 
gram in this work is to investigate the effect of positional 
disorder itself on transport characteristics. Although the 
treatment is nonperturbative with no local order in the 
locations of the sites in the amorphous systems we ex- 
amine, the decay length scale I of the hopping integral 
serves to parameterize the disorder strength. For large 
I (small 7), the tunneling is long-ranged and connects 
sites to many neighbors, effectively averaging over the 
hopping rates to and from many orbitals and thereby 
somewhat muting the effect of positional disorder On the 
other hand, if / is small relative to the typical inter-site 
spacing p~ x l D (i.e. for large 7), the tunneling is prepon- 
derantly to the very closest neighbors; moreover, even 
small fluctuations in the locations of nearest neighbors 
impact charge transport through the site to a significant 
degree via the exponential dependence of the hopping 
integral. Hence, disorder is in a sense amplified as 7 is 
increased and muted when 7 is small. 

In calculating the tight binding wave functions, we ex- 
amine a L D supercell, where many system sizes are con- 
sidered in order to perform finite size scaling and de- 
termine the degree to which eigenstates are localized in 
the bulk limit. Energies obtained by diagonalizing the 
tight binding (Hermitian) Hamiltonian are used to con- 
struct the global density of states (DOS), while the eigen- 
states themselves are retained for analysis to character- 
ize the electronic wave functions with respect to local- 
ization with the aid of a single quantity known as the 
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Inverse Participation Ratio (IPR), ^]|V'i| 4 /( 5Zl^l 2 ) ■ 
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The participation ratio shows distinct behavior depend- 
ing on whether the wave function is confined to a small 
volume or spread out over a larger region, and hence 
more extended in character. While in the former case 
the IPR is finite and tends to be relatively large, the par- 
ticipation ratio is smaller for broader electronic states, 
and approaches zero as the wave function becomes spread 



over a bulk system in the case of a genuinely extended 
state. This dichotomy for extended vis a vis localized 
states makes the IPR a useful diagnostic parameter in 
the context of theoretical calculations where often the 
participation ratio provides information as to the extent 
of localization of carrier wave functions in specific lo- 
cations in an amorphous geometry or in certain energy 
ranges of a band structure^— . 

In this work, we aim to determine for very large ID, 
2D, and 3D amorphous systems the prevalence of ex- 
tended states, and we construct phase portraits showing 
where wave functions are localized and where states are 
extended. The random character of the disorder pre- 
cludes a direct observation of the evolution of the char- 
acteristics of electronic states as the system size L is in- 
creased. However, there remains the possibility of deter- 
mining in a statistical sense the localization characteris- 
tics of the amorphous medium by calculating the Inverse 
Participation Ratio histogram. Shifts in the weight of 
the IPR probability distribution with increasing L pro- 
vide information as to how many of the electronic states 
ultimately are localized and what portion are genuinely 
extended. For a specific hopping integral decay parame- 
ter 7 we find the status of electronic states with respect to 
localization to be determined exclusively by the energy 
eigenvalue, with no situation arising in which localized 
and extended states exist simultaneously for the same 
infinitesimal energy interval. With a finite size scaling 
analysis of the Participation ratio histogram calculated 
for various system sizes, we extrapolate to the bulk limit 
and determine in a rigorous fashion energies where elec- 
tronic states are localized, and eigenenergies supporting 
extended states. 

By repeating the calculation for different 7, we obtain 
phase portraits showing regions of localized and extended 
wave functions. The IPR probability distributions rep- 
resent an intermediate stage in the determination of en- 
ergy ranges corresponding to localized or extended states 
in the bulk limit, but are of interest in themselves and 
highlight salient qualitative trends. A participation ra- 
tio histogram which does not change either in its shape 
or position with increasing system size signifies that all 
states encompassed in the distribution are localized, and 
the IPR distribution may be regarded as a bulk charac- 
teristic where L ;§> £ with £ being the localization length 
scale of the eigenstates. At the opposite extreme is a 
Participation ratio histogram which shifts in its entirety 
toward smaller IPR values with increasing L. The steady 
transfer of a large share of statistical weight toward even 
smaller participation ratios implies the electronic states 
are preponderantly extended. A possibility interpolat- 
ing between the two extremes arises if a part of the IPR 
distribution converges as would be expected for a set of 
localized states while at the same time a portion of the 
histogram separates from the main envelope of the dis- 
tribution and is conveyed toward lower IPR values. 

To sample realizations of disorder from the appropri- 
ate statistical distribution, it is important to generate 
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random configurations of atoms in an unbiased way. For 
stochastic input, we use a Mersenne Twister algorithm 
to minimize correlations among successively generated 
random numbers and to ensure the period of the pseu- 
dorandom sequence far exceeds the quantity of random 
numbers used over the course of the simulations. For con- 
tinuously distributed hopping sites, the number N of sites 
in the simulation volume must in general vary from one 
sample to the next due to statistical fluctuations. The in- 
teger value closest to the mean occupancy N av = pL D is 
a convenient initial choice, and random variations in the 
number of particles in the simulation volume are taken 
into account with a sequence of stochastically driven at- 
tempts either to raise or lower N . The latter are part of 
an importance sampling scheme similar to that used to 
derive the Metropolis Criterion^ at the heart of Monte 
Carlo simulations which sample the Boltzmann distribu- 
tion in the calculation of thermodynamic variables. 

To determine the probability of N sites in v = L D , 
we divide v into M sub-volumes of equal size where 
Av = v/M. For large M the likelihood of multiple occu- 
pancy in any of the sub- volumes is very small relative to 
the chance of having one or zero sites in a subdivision; 
in the small Av limit, the single occupancy probability 
is pAv, with (1 — pAv) being the complementary likeli- 
hood of null occupancy. Hence, the probability the en- 
tire system is devoid of hopping sites is (1 — pv/M) M , 
which becomes e~ pv for M — > oo. For single occupancy, 
adopting a prefactor M to take into consideration that 
the site may reside in any of the M sub- volumes, yields 
M(pv/M)(l - pv/M) M ~ r , which becomes pve pv in the 
6v — > limit. Similar logic is used for general case, and 
the probability for having exactly N sites in the simu- 
lation volume is P(N) — e~ pv (pv) /N\ where AM is a 
combinatorial factor to compensate for multiple count- 
ing. 

To generate a realization of disorder, a succession of 
attempts (a number of moves in the vicinity of N av is 
sufficient to achieve ergodicity) is made to raise or lower 
the occupancy number N, where the choice to increase 
or decrease N is randomly determined. For increments 
from N to N + 1, the relevant criterion is the probabil- 
ity ratio r + = p(N + l)/p(N) = pv/(N + 1), and the 
change is accepted is X r < r+, where X r is a random 
number sampled uniformly from the interval [0, 1]. Sim- 
ilarly, decreasing N to N — 1 occurs if X r < r_ where 
r_ = p(N — l)/p(N) = N/ pv. With N properly sampled, 
D Cartesian coordinates for each site location are chosen 
independently (and at random with uniform probability 
density) from the interval [0,L]. 

The coordinates for each of the N sites enter in the 
construction of the Hamiltonian matrix (in the context 
of the tight binding model given in Eq. [I}, which is di- 
agonalized for the eigenenergies and eigenstates; periodic 
boundary conditions are implemented to mitigate finite 
size effects. For the purpose of the DOS calculations, 
5 x 10 5 energy eigenvalues are sampled; in obtaining the 
IPR statistical distributions, where eigenstates are used 



to calculate participation ratios with a concomitant in- 
crease in the computational burden, 10 5 wave functions 
are retained. 

In each section of this work, we present and discuss 
results for ID, 2D, and 3D systems; for each dimension- 
ality, we examine a range of tunneling matrix element 
decay parameters 7 which in a sense govern the strength 
of the disorder. In Section II, we examine the Density of 
States for energy eigenvalues. Inverse Participation Ratio 
(IPR) statistical distributions are discussed in Section III 
with IPR histograms displayed for various system sizes. 
In Section IV, the channel averaged participation ratios 
are shown; we argue that IPR channel averages are rep- 
resentative of a particular energy and hence may be used 
to determine how the IPR scales with system size. Ul- 
timately, in Section V, a finite size scaling analysis of 
the participation ratio channel averages is used to ex- 
trapolate to the thermodynamic limit and thereby yield 
the mean bulk IPR. The latter, used to construct phase 
portraits showing where states are extended or localized 
for different values of E and 7, indicates for D = 3 a 
zone of extended states in the midst of a region of lo- 
calized states and flanked by two mobility boundaries. 
Although a finite fraction of the states are extended even 
for large 7 where disorder fluctuations are very impor- 
tant, the width of the interval of energies 111(7) where 
states are extended decreases rapidly for large to mod- 
erate 7, with an asymptotically exponential dependence 
1/7(7) ~ e~ Al . Conclusions are presented in Section VI. 



II. ENERGY DENSITY OF STATES 

The DOS, which in principal is a continuous func- 
tion f(E) of energy for an amorphous material, may 
only be rendered to an approximate degree in a finite 
calculation. Accordingly, we partition the energy inter- 
val between the ground state and the uppermost excited 
state into a finite though reasonably large number of sub- 
intervals or "bins". The augmented resolution achieved 
with an increase in the total number of partitions is coun- 
terbalanced with a rise in the magnitude of statistical 
fluctuations. To strike a suitable balance among detail 
and noise, 500 divisions are used in preparing DOS his- 
tograms; with a total of 5 x 10 5 eigenvalues for each DOS 
curve, choosing 500 bins still allows on average for 1,000 
data point for each partition. 

We examine DOS curves for ID, 2D, and 3D systems, 
for which a salient common characteristic is a rapid con- 
vergence of the energy eigenvalue statistical distributions 
with respect to the system size L. In fact, as may be 
seen for representative cases in Fig. [TJ the DOS curves 
overlap very closely and approach the bulk limit as long 
as L > maxfp" 1 / 15 ,^ with the hopping range I and the 
typical interparticle separation p~ x l D (unity in our treat- 
ment) being relevant length scales. 

Another important trend seen for all dimensionalities 
under consideration is a systematic shift of the DOS sta- 
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FIG. 1: (Color Online) Density of states plotted for (top to 
bottom) D = 1, D = 2, and D — 3, and for large to small 
7 (left to right). DOS curves are plotted for various systems 
sizes L. 
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FIG. 2: (Color Online) Density of States curves are shown 
together for a range of decay rates 7, with, clockwise from 
left, D = 1, D = 2, and D = 3. 
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FIG. 3: (Color Online) Inverse Participation Ratio profiles 
graphed for various system sizes L for hopping integral decay 
rates ranging from 7 = 3.0 in panel (a) to 7 = 0.5 in panel 
(d) for ID systems. 



tistical weight toward lower (i.e. more negative) energies 
with decreasing 7, or increasing tunneling matrix element 
range I — 7 -1 shown. In Fig.[2]DOS curves are displayed 
for D = {1, 2, 3} for a broad range of 7 values with re- 
sults (calculated for large L) converged with respect to 
the size of the system. 

Although a shift of statistical weight toward negative 
energy values occurs for each dimensionality, the evolu- 
tion of the DOS curves for ID systems differs from that 
seen for D = 2 and D — 3. For D = 1, the energy proba- 
bility distribution has a single maximum centered about 
E = for 7 > 1. However, with increasing / = 7 _1 , a 
secondary peak appears at a negative energy, gaining am- 
plitude at the expense of the probability density near the 
E = maximum; for large enough /, the latter vanishes 
with only the negative energy peak remaining. 

For two and three dimensional systems, probability 
density is transferred from E = to negative energies 
as in the case D = 1. However, instead of an interme- 
diate transition to a dual peak profile, the DOS curves 
remain unimodal as the single maximum migrates con- 
tinuously toward lower energies, arriving at a negative 
energy for 7 <C 1. 



III. PARTICIPATION RATIO STATISTICAL 
DISTRIBUTIONS 

The Inverse Participation Ratio (IPR) is a compact 
single parameter gauging the degree to which electronic 
states are localized or extended, with the IPR tending 
to zero for L — > 00 for bulk extended states. While our 
aim in this work is to extrapolate to the bulk limit in 
a quantitative fashion, information may also be gleaned 
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FIG. 4: (Color Online) Inverse Participation Ratio profiles 
graphed for various system sizes L for hopping integral decay 
rates decay rates ranging from 7 = 3.0 in panel (a) to 7 = 0.5 
in panel (f) for 2D systems. 
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FIG. 5: (Color Online) Inverse Participation Ratio profiles 
graphed for various system size L for hopping integral decay 
rates 7 ranging from 7 = 3.5 in panel (a) to 7 = 1.0 in panel 
(f) for 3D systems. 



at a qualitative level when participation ratio curves are 
juxtaposed for a range of system sizes. Since the IPR 
may vary by several orders of magnitude over the full 
gamut of system sizes L under consideration, it is often 
more prudent to exhibit log 10 (IPR) in lieu of the raw 
participation ratios. 

The IPR histograms are created by dividing the inter- 
val along the log 10 (IPR) abscissa into a suitable number 
of bins. With the availability of 10 5 participation ratios, 
the choice of 200 partitions provides a reasonable mea- 
sure of resolution while keeping statistical fluctuations at 
reasonable levels. IPR distributions are shown in Fig. |31 
Fig. |H and Fig. [S] for D = 1, D = 2, and D = 3 respec- 
tively. 

In the ID case, the system sizes progress geometrically, 
doubling from L — 100 to L = 6400. With successive 
doublings of L, the IPR histograms invariably converge 
and cease to evolve with increases in the system size, a 
characteristic consistent with the localization of all wave 
functions in the thermodynamic limit. When converged 
in L to an IPR distribution appropriate to the bulk limit, 
histograms are dominated by cusp-like peaks near the 
upper limit of the IPR range. The latter characteristic 
is a hallmark particular to ID systems and is evident 
whether the hopping integral length scale / is large or 



small relative to the separation between sites. 

The interpretation from the D = 1 participation ra- 
tio distributions that all wave functions are localized is 
in a sense not surprising. With the theoretical frame- 
work of Anderson localization having been introduced 
more than 50 years ago^, a significant body of work 
(both in experiment and theoretical calculations) has ex- 
amined the tendency for random potentials to localize 
electronic states very effectively in one dimensional sys- 
tems, even for weak random potentials. Moreover, the 
availability of cold atom traps with coherent quantum 
states where the underlying one dimensional potential 
may be tailored in a variety of ways has made possi- 
ble the study of localization properties of ID systems in 
a controlled manner. In this vein, a direct experimen- 
tal observation of localization has recently been achieved 
in a Bose-Einstein condensate with the random (diag- 
onal) potential set up by a laser speckle^ with results 
in accord with theoretical descriptions- 31 . Bichromatic 
aperiodic potentials are not purely random uncorrelated 
disorder, but nonetheless have been found in experiment 
and theoretical analysis^— to be very effective in local- 
izing quantum states. In this work, we show using finite 
size scaling analysis that off-diagonal disorder inherent in 
one dimensional amorphous systems leads to the localiza- 
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tion of all electronic states in the absence of random site 
potentials. 

For D — 2, the evolution of the participation ratio 
with increasing L depends on the hopping integral decay 
length I. IPR histograms for various 7 values are shown 
in the six panels of Fig. 2] as in the ID case, successive 
values of L differ by a factor of 2, and participation ratio 
distributions for L = {5,10,20,40,80} appear together 
on the same graph. For short hopping ranges I = 7 _1 , 
there are maxima in the high IPR regime which sub- 
sume a large share of the total statistical weight; apart 
from small vacillations, the distributions have a unimodal 
envelope, as seen in IPR histograms calculated for one 
dimensional systems. Moreover, the convergence of his- 
tograms with respect to L is an indication of the local- 
ization of all states in the thermodynamic limit. 

For longer decay lengths I (i.e. for 7 < 1.5) partic- 
ipation ratio distributions obtained for 2D amorphous 
systems lack the unimodal sharply peaked structure seen 
for small I and unlike the I <C 1 counterparts, seem not 
to converge with respect to increases in L. Instead, his- 
tograms become bimodal as statistical weight is trans- 
ferred to the left, toward the lower IPR regime. 

The migration of statistical weight toward lower IPR 
values and more extended character, is indicative of the 
possible existence of extended states, and the effect is 
particular striking for the smallest decay rates 7 = 0.75 
and 7 = 0.5. The systematic transfer with increasing L 
of probability density to smaller IPR values is manifest 
as a steady leftward shift of the trailing (low IPR) edge 
and of the histogram with the nearby maximum also born 
leftward. In fact, the shift both of the peak and the left- 
most front is constant in magnitude each time the system 
size L is doubled; since the abscissa is log 10 (IPR), the 
latter trend corresponds to a power law scaling for 
the left trailing edge and low IPR peak. 

Participation ratio histograms for 3D systems are dis- 
played in Fig. [5] As in the large / limit for 2D systems, 
a portion of the participation ratio statistical weight 
shifts systematically to lower IPR values. The size of 
the emerging peak and the amount of probability den- 
sity which migrates leftward increases with /, even as the 
cusp-shaped maximum in the low IPR regime decreases 
in amplitude and overall statistical weight. 

A noteworthy feature for I < 1 is the robustness of 
the probability density contained in the leftward shift- 
ing peak. That the packet moves systematically toward 
lower IPR values without leaving behind any statistical 
weight is compatible with the existence of a finite frac- 
tion of extended states in the bulk limit. Thus, for 3D 
amorphous systems, there is a dichotomy in the way the 
distribution changes, where the leftward shifting proba- 
bility density corresponding to extended wave functions 
contrasts with the localized states encompassed in the 
high IPR region of the histogram, which ceases to evolve 
with increasing L. With increasing tunneling matrix ele- 
ment range I, the balance shifts in favor of the extended 
states as more and more statistical weight is swept into 
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FIG. 6: (Color Online) Log-Log graphs of channel averaged 
IPR versus L for ID systems in the case of a 100 channel 
scheme. Panel (a) and panel (b) correspond to hopping inte- 
gral decay rates 7 = 1.0 and 7 = 0.5 respectively. 
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FIG. 7: (Color Online) Log-Log graphs of channel averaged 
IPR versus L for 2D systems in the case of a 100 channel 
scheme with hopping integral decay rates ranging from 7 = 
2.5 in panel (a) to 7 = 0.5 in panel (f). 



the peak moving toward lower IPR values. Nonetheless, 
as embodied in the part of the distribution which does 
not change with increasing L, a finite fraction of the wave 
functions are localized even for very small values of the 
decay rate 7. 
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FIG. 8: (Color Online) Log-Log graphs of channel averaged 
IPR versus L for 3D systems in the case of a 100 channel 
scheme with hopping integral decay rates ranging from 7 = 
4.5 in panel (a) to 7 = 0.75 in panel (f). 



IV. THE CHANNEL AVERAGED 
PARTICIPATION RATIO 

With long-range positional order absent in the amor- 
phous systems we examine, the eigenstate energy is the 
only good quantum number available, and in this work 
we show that the energy eigenvalue is a unique determi- 
nant as to whether electronic states are either extended 
or localized for a particular energy with the simultaneous 
presence of localized and extended states ruled out as a 
possibility. 

Determining if a wave function tp is localized or ex- 
tended entails calculating participation ratio statistics for 
many L values and using finite size scaling to access the 
L — > 00 limit, with a vanishing IPR in the bulk limit a 
hallmark of extended states. The random character of 
the disorder precludes the study of the evolution of indi- 
vidual states with increasing system size, and instead we 
must analyze aggregates of wave functions across a range 
of system sizes. Electronic states are parameterized by 
energy eigenvalues, and one possible choice is to partition 
the states into channels centered about uniformly spaced 
energies E which encompass an energy interval SE nar- 
row enough to capture information specific to wave func- 
tions with energies very close to E , but broad enough to 
suppress statistical fluctuations. 



The DOS statistical distributions shown in Fig. [2] are 
sharply peaked, with a rapid decrease in the probability 
density away from the maximum. In a practical sense, 
the DOS heterogeneity poses a challenge for a scheme 
where the global energy range AE between the ground 
state energy E m i n and the highest excited state energy 
£ max is partitioned into small intervals SE of uniform 
size; statistical fluctuations will plague channels far from 
the DOS maximum where the statistical weight for eigen- 
states is sharply reduced. 

In lieu of energy, to circumvent the problem of non- 
uniform statistics, states are labeled with the normalized 
energy eigenvalue rank f. For finite systems, the rank 
number is assigned by calculating energy eigenvalues and 
corresponding wave functions for a large number of con- 
figurations of disorder. The normalized rank number is 
f = r/N where N is the total number of states and r 
is the global eigenvalue rank within the large aggregate. 
We find n = 100 yields channel widths Sf sufficiently 
narrow for channel averages to be representative of the 
normalized rank f at the center of the channel, yet broad 
enough to provide sufficient statistics for analysis; a par- 
allel calculation for n — 50 yields results in quantitative 
agreement with the n = 100 scheme, direct verification 
that n — 100 is large enough to avoid systematic admix- 
ture effects from the finite channel width Sf. 

Channel averaged IPR results appear in Fig. [51 Fig. [71 
and Fig. [5] for ID, 2D, and 3D systems, where the hor- 
izontal axis is log 10 (L) for each case. The participation 
ratio results for the one dimensional systems are in a 
sense most readily interpreted. IPR curves initially de- 
crease for small system sizes, but quickly level out and 
approach asymptotically finite values even for the small 
decay constant 7 = 0.5. The latter phenomenon is con- 
sistent with the convergence of the global IPR probabil- 
ity distribution with increasing L, interpreted as a sign 
that all wave functions are localized in one dimensional 
systems irrespective of the hopping decay parameter 7. 

For D = 2, there is a bifurcation in the way the IPR 
channel averages vary with increasing L with the precise 
behavior determined by the hopping range I — 7 _1 . For 
large 7 (i.e. especially for 7 = 2.0 and 7 = 2.5), the IPR 
traces level out and approach finite asymptotic bulk val- 
ues corresponding to localized states, much as occurs for 
ID amorphous systems. On the other hand, in the case of 
smaller 7, where there is a greater likelihood that states 
with extended character are supported, it is not defini- 
tively conclusive that the participation ratio curves ap- 
proach asymptotically finite values for very large system 
sizes. Broadly speaking, for D = 2 IPR traces are most 
likely to level out and tend to a finite value for smaller 
f, or for energies near the ground state. For 7 < 1.0, 
there are channel average curves which in principal may 
keep a finite negative slope as L — ¥ 00, which on a log- 
arithmic scale is tantamount to a vanishing IPR in the 
thermodynamic limit. 

Notwithstanding the persistently downward slope, an 
unambiguous determination that the monotonically de- 
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FIG. 9: (Color Online) Intra-channel standard deviations of 
log 10 (IPR) plotted for the D = 1 case. Graphs in panels (a) 
and (c) correspond to relatively large 7 = 1.5, while plots in 
panels (b) and (d) are calculated for a more gradual hopping 
integral decay, 7 = 0.5. Graphs on the left show a versus 
log 10 (L) for various normalized channel numbers f, while the 
rightmost panels are plots of <r[log 10 (IPR] with respect to the 
rank r for large L. 



FIG. 10: (Color Online) Intra-channel standard deviations of 
log 10 (IPR) plotted for the D — 2 case. Graphs in panels (a) 
and (c) correspond to relatively large 7 = 2.0, while plots in 
panels (b) and (d) are calculated for a more gradual hopping 
integral decay, 7 = 0.5. Graphs on the left show a versus 
log 10 (L) for various normalized channel numbers f, while the 
rightmost panels are plots of cr[log 10 (IPR] with respect to the 
rank f for large L. 



creasing channel averages represent extended states is 
hampered by the upward concavity in most of the IPR 
curves, where except for the very highest f (i.e. for 
f = 0.95 where the concavity appears to be neutral), 
there is a progressive reduction in the magnitude of the 
downward slope which could eventually cause the chan- 
nel averaged participation ratio to level out at a finite 
value. To determine in a rigorous way if the participa- 
tion ratios vanish for L — > 00 or instead approach a finite 
value, a finite size scaling analysis, described in Section 
V, is needed. 

As in the case of the 2D systems, for D = 3, there are 
in a broad sense two ways in which IPR curves scale with 
L. For very large decay rates (e.g. 7 = 4.5), the partici- 
pation ratio channel averages seem to level out and tend 
to finite values. On the other hand for smaller 7, instead 
of the upward concavity seen in the case of two dimen- 
sional systems, many of the curves are concave downward 
with the rate of decrease of the channel averages increas- 
ing with L. With the downward slope becoming greater 
instead of showing signs of faltering as in D = 2, it is 
possible to conclude without further analysis that the 
participation ratio tends to zero in the bulk limit. 

With channel averages only providing the mean partic- 
ipation ratio, it is important to be certain the IPR values 
obtained in this manner represent the characteristics of 
all of the states encompassed in a channel. The most 



vital task in this vein is to rule out the possibility of the 
coexistence of localized and extended states for a specific 
channel index r (or energy E) . To tackle this question, we 
calculate the standard deviation a of log 10 (IPR) across a 
channel with results appearing Fig. |H] for D = 1, Fig. [TU1 
for D = 2, and Fig. [TT] for D — 3. With increasing L, the 
quantity cr[log 10 (IPR)] will either diverge or tend to a fi- 
nite value. While a divergence indicates a dual presence 
of extended and localized states (coexistence) , a standard 
deviation of participation ratio logarithms which tends to 
a finite value for L — > 00 is an indication that there are 
either localized or extended wave functions, but not both 
(exclusivity) . 

If only localized states are present, the the standard 
deviation will converge to a limiting value (cr is in gen- 
eral finite even in the thermodynamic limit since some 
wave functions are more spread out than others due to 
disorder fluctuations) and become a bulk characteristic 
when L ^> £, where £ is the largest decay length scale 
of the broadest states. On the other hand, for a suite 
of purely extended states, the IPR for each wave func- 
tion will tend to zero with a concomitant divergence of 
log 10 (IPR). Although the magnitude of the ensemble 
average (log 10 (IPR)), which is negative in sign, grows 
without bound for increasing L for a rank f or energy 
E supporting extended states, it is nevertheless not clear 
that a diverges since the standard deviation measures the 
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FIG. 11: (Color Online) Intra-channel standard deviations of 
log 10 (IPR) plotted for the D — 3 case. Graphs in panels (a) 
and (c) correspond to relatively large 7 = 3.0, while plots in 
panels (b) and (d) are calculated for a more gradual hopping 
integral decay, 7 = 1.0. Graphs on the left show a versus 
log 10 (L) for various normalized channel numbers f, while the 
rightmost panels are plots of <r[log 10 (IPR] with respect to the 
rank f for large L. 



spread in the participation ratios, not their magnitude. 
In the global IPR histograms shown in Fig. U and Fig. [5] 
for D = 2 and D = 3, the fact that packets of probability 
density which migrate lower log 10 (IPR) with increasing 
L tend to maintain their shape and profile even as the 
participation ratio decreases by several orders of magni- 
tude suggests a may remain finite even as the average 
(log 10 (IPR)) diverges. 

In the coexistence scenario, with a simultaneous pres- 
ence of extended and localized wave functions, there 
will be a divergence in cr[log 10 (IPR)] since the packet 
of probability density for localized states remains fixed 
with increasing L while the peak corresponding to states 
with extended character is conveyed toward more nega- 
tive log 10 (IPR). Since the standard deviation provides 
a measure of the increasing separation of the two peaks, 
cr[log 10 (IPR)] must diverge as L — > 00 if localized and ex- 
tended electronic states are present for the same energy 
E or rank f. 

To determine if for any combination of system pa- 
rameters a converges to a finite value (exclusivity) 
or diverges (coexistence), we plot cr[log 10 (IPR)] versus 
log 10 (£) in the left panels of the graphs in Figs.li-JlFigs.ITUl 
and Figs. [11] In addition, the right panels show a 
with respect to f for large L to provide a glimpse of 
cr[log 10 (IPR)] for large L. The er curves for the cases 
D = 1, and D = 2, and D = 3 are concave downward 
in the large L regime, and seem to level out and tend 



to finite values. In addition, for the largest system size 
considered in this work and plotted in the right panels 
of Fig. U Fig. [TUl and Fig. QTJ <r[log 10 (IPR)] are low in 
magnitude, and do not exceed cr[log 10 (IPR)] = 0.35 for 
which the spread of the participation ratios is at most a 
factor of two. 



V. PARTICIPATION RATIOS IN THE 
THERMODYNAMIC LIMIT 

The primary aim in the finite size scaling analysis is to 
calculate IPRo, the channel averaged participation ratio 
in the thermodynamic limit. For the IPR dependence 
for moderate to large L, we use a power law formula, 
IPR(i) = IPR + ot\L-P + a 2 L~ s where IPR is the 
participation ratio in the bulk limit, j3 is the leading or- 
der scaling exponent, and 5 is the exponent for the next 
to leading order contribution to scaling; a\ and ai are 
amplitudes. We obtain the parameters IPRo, P, S, 011, 
and «2 in a nonlinear least squares calculation by min- 
imizing the sum of the square of the relative differences 

A — 1 rv^ m /IPR CA (L;)-IPR LSF (L,) N2U/2 /• -4.1, 
A LSF = £Ei=l( IPR^A (£i) T ] ' (!- e - Wlth 

L\ and L m the smallest and largest systems examined, 
respectively) of the data gleaned from the IPR channel 
averages IPR CA (L) and the theoretical scaling expression 
IPR LSF (L). To find the optimal fit, we use a stochastic 
algorithm with the quantity Alsf treated as an "energy" 
to be minimized by randomly perturbing IPRo, ampli- 
tudes ai and c*2, and exponents /3 and S. Only Monte 
Carlo moves which decrease Alsf (and hence incremen- 
tally improve the fit) are accepted, and the IPR(L) pa- 
rameters are suitably converged after 4 x 10 5 attempts to 
shift the five unknown parameters in IPR(L) in a stochas- 
tic fashion. 

Participation ratios extrapolated to the bulk limit are 
displayed for D = 1, D = 2, and D = 3 in Fig. QJJ The 
main graphs are calculated for 100 channel partitions, 
whereas the insets show results for a 50 channel scheme. 
The good agreement among the 100 and 50 channel bulk 
IPR values indicates convergence of the bulk IPR values 
with respect to the number of channels. 

Bulk limit Participation ratios for ID systems are plot- 
ted in panel (a) of Fig. [T2] Notwithstanding a nonmono- 
tonic variation of IPR in the normalized channel number 
f the D = 1 results are finite in all cases, dipping only 
slightly below IPRo =0.1 even for hopping integral de- 
cay rates as low as 7 = 0.5, and thus are indicative of 
localization of all wave functions irrespective of the size 
of the hopping range I. 

Extrapolated IPR results for D = 2 shown in panel 
(b) of Fig. [T2] are non-monotonic, generally decreasing 
precipitously near r = and ultimately recovering in the 
vicinity of r = 1. With increasing I, the participation 
ratio trough becomes broader and deeper, while the re- 
covery to higher IPR values is muted. Eventually, for 
7 < 1, the minimum reaches the horizontal axis, where 
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IPRo = 0, and there is no return to higher participation 
ratios near f = 1; hence the f axis is effectively divided 
into two regions where wave functions are localized below 
a threshold r c and extended above it. 

The D = 3 bulk IPR curves in panel (c) of Fig.QJdiffcr 
in significant ways from the corresponding participation 
ratios calculated either for ID or 2D systems. In general, 
looking from r = to f = 1, the bulk IPR initially makes 
a sharp descent until intersecting the horizontal axis. The 
slope downward decreases with the hopping range /, and 
the normalized channel number where the bulk IPR be- 
comes zero retreats to smaller f with increasing scale 
For 7 < 1, there is only a single mobility boundary; the 
IPR falls to zero, coincides with the abscissa, and does 
not climb again to finite values. 

On the other hand, although for larger 7 the extrap- 
olated participation ratio curves initially descend as in 
7 < 1 and eventually become zero, in the upper f range, 
the IPR begins to increase and is again finite. Thus, for 
7 > 1, two mobility boundaries separate an intermediate 
range of energies where the participation ratio vanishes 
(corresponding to extended states) from energy regimes 
where wave functions are localized. The width of the in- 
terval where the IPRo = constricts width increasing 7, 
though even for quite large values of the decay constant 
(i.e. for 7 = 5.5), the extrapolated IPR seems at least 
to briefly touch the horizontal axis before rising again to 
finite values with a locally "v" shaped profile. Participa- 
tion ratio curves corresponding to smaller 7 values have 
broader minima, but also have local symmetry about the 
same point in the vicinity of f — 0.58. 

We use bulk participation ratios to construct phase di- 
agrams for the tight binding wave functions; for the ID 
case, there is only a single phase with the electronic states 
are localized in all cases, while phase portraits for D = 2 
(shown in Fig. [13] and Fig. [T4l) and D — 3 (displayed in 
Fig. US and Fig. (THl) show areas of extended states and 
phases in which all of the wave functions are localized. 
For the 2D case, the wedge of extended states, broadest 
for 7 = 0.5, tapers with increasing 7 and vanishes al- 
together in the vicinity of 7 — 1.5, as indicated by the 
broken lines which are extrapolations. 

For D = 3, the swath of extended states is quite broad 
for 7 ~ 1.0, encompassing most of the eigenstates and 
the associated range of eigenstates energies. As in the 
2D smaller portion of the wave functions are ex- 

tended with increasing 7. for moderate to large hopping 
integral decay lengths (i.e. for 7 > 3.0), the band of ex- 
tended states rapidly constricts with the extended phases 
asymptotically symmetric about f = 0.58. Fig. [T5]is the 
phase diagram with the vertical axis rendered in terms 
of energies instead of normalized rank r. Again, the in- 
terval of extended states becomes sharply narrower with 
increasing 7. 

Although the extended state phase persists for large 
values of 7 (i.e. even for 7 > 5.0), the decrease of the 
width is very rapid, and we examine the possibility that 
the decrease may be exponential in the hopping integral 
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FIG. 12: (Color Online) Inverse Participation Ratio profiles 
graphed for strong decay rates 7 [ranging from 7 = 3.5 in 
panel (a) to 7 = 2.0 in panel (d)] for ID systems. 
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FIG. 13: (Color Online) Phase portrait for 2D systems with 
the hopping decay parameter 7 on the abscissa and the chan- 
nel number f on the vertical axis. The broken lines are ex- 
trapolated parts of the phase boundary. 
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FIG. 14: (Color Online) Phase diagram for the D = 2 case 
with the hopping decay parameter 7 on the abscissa and the 
eigenstate energy (in electron volts) on the vertical axis. The 
broken lines indicates extrapolated parts of the phase bound- 
ary. 
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FIG. 15: (Color Online) Phase portrait for 3D systems with 
the hopping decay parameter 7 on the abscissa and the chan- 
nel number f on the vertical axis. 
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decay rate 7. Fig. [T7] displays a graph of base ten log- 
arithm of the width of the phase where eigenstates are 
extended; the main plot shows log 10 [w(f)] with respect 
to 7, while the inset is a plot of \og 1Q [w(E)} versus 7. The 
asymptotically linear dependence of the logarithm of w, 
which is seen whether one considers the the width w(E) 
of the energy interval or w(f) of the normalized rank, is 
consistent with an exponential dependence w oc e~ Al for 
moderate to large 7. 



VI. CONCLUSIONS AND FUTURE WORK 

We have calculated the energy Density of States and, 
using the Inverse Participation Ratio, we have examined 
the characteristics of electronic states in amorphous sys- 
tems in one, two, and three dimensions for hopping ma- 
trix elements which decay exponentially in the separation 



FIG. 16: (Color Online) Phase diagram for the D = 3 case 
with the hopping decay parameter 7 on the abscissa and the 
eigenstate energy (in electron volts) on the vertical axis. The 
graph inset is a closer view of the bifurcation region. 



distance between neighboring sites in the context of a 
tight binding model. We have calculated global IPR sta- 
tistical distributions, which have a rich multimodal struc- 
ture for systems in two and three dimensions in contrast 
to the simple sharply peaked profiles consistently seen for 
D = l. 

Partitioning wave functions according to the normal- 
ized eigenvalue rank f has yielded channel averaged par- 
ticipation ratios, which we have shown to be represen- 
tative of wave functions for a specific energy E or nor- 
malized rank f. By applying finite size scaling to each 
of the channel averages, we have obtained participation 
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FIG. 17: (Color Online) Logarithm of the width of the ex- 
tended state region plotted versus the decay parameter 7. 
The main graph shows the normalized rank width w(r), and 
the inset plot displays the width of the corresponding energy 
interval w(E). 

ratios in the thermodynamic limit; using the latter, we 
have constructed phase portraits of the eigenstates with 
respect to localization. 

In ID, wave functions are strongly localized in all cases, 
whereas results for 2D indicate the presence of a critical 
decay parameter j c , with localization for 7 > 7 C and 
existence of extended states below j c . In the 3D case, 
extended states also are admitted, and the wave func- 
tions with extended character occur even for quite large 
7 values, although the interval of energies supporting ex- 
tended states diminishes with 7, asymptotically scaling 
as e~ Al with increasing 7. The swath of extended states 
is flanked by regions where wave functions are localized, 
with the two interfaces interpreted as mobility bound- 



aries. 

In future studies, disorder schemes will be considered 
in which the severity of the disorder is tuned from mild 
to quite strong by perturbing a regular periodic crystal 
lattice and introducing random perturbations S in the po- 
sitions of the hopping sites. The disordering shifts S may 
be introduced, e.g., from a Gaussian distribution with 
a RMS magnitude a. Among the salient germane ques- 
tions to be investigated in this manner is whether there is 
a threshold in typical displacement magnitudes where ex- 
tended states may survive in D — 1 and D — 2 if random 
displacements in atomic positions are sufficiently small in 
relation to the crystal lattice constant. Given the fragility 
of extended character in D = 1 , one might predict that 
even a small random perturbation in the site positions 
from a periodic configuration might induce localization 
in a one dimensional lattice. On the other hand, in two 
dimensions it may be that there is a perturbation level 
beyond which the disordering influence causes most or 
all of the states to be localized, with predominantly ex- 
tended character below the perturbation threshold. 

In the present study, we have concentrated on short- 
range couplings, as might be appropriate in an exchange 
type coupling scheme. Nonetheless, it would be useful 
also to examine a power law decay to see if the severity 
of localization effects are reduced in one dimension, and 
if bona fide extended states exist under these conditions. 
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